home *** CD-ROM | disk | FTP | other *** search
/ Power Programmierung / Power-Programmierung (Tewi)(1994).iso / magazine / aijournl / 1988_08 / kohonen.c < prev    next >
Text File  |  1988-07-22  |  16KB  |  455 lines
  1. /*  Kohonen Feature Map simulation */
  2. /*********************************************************
  3.   Kohonen Network
  4.   Written by Maureen Caudill 
  5.   in Lightspeed C on a Macintosh
  6.   
  7.  See the August, 1988 AI Expert, Neural Network Primer, Part 4 
  8.   for a discussion of this network.
  9.   
  10.   This program simulates a Kohonen feature map layer.  A total of NUMITERS 
  11.   2-dimensional data patterns are presented to the network.  After each
  12.   SAVECOUNT patterns, a data file is written which saves the current 
  13.   state of the weight vectors for the network.  These snapshots, when 
  14.   properly plotted, will generate the topological map of the system at each
  15.   point in time.  See Kohonen Plot.C for source code detailing how to draw
  16.   the feature map.  
  17.   
  18.   This program is generic to any computer system.  Kohonen Plot.C 
  19.   uses Macintosh-specific code to draw the graphics.  If another computer
  20.   is used, individual graphics routines will have to be written by the
  21.   user to generate on-screen graphics.
  22.   
  23.   ⌐Adaptics
  24.   16776 Bernardo Center Drive, Suite 110B
  25.   San Diego, CA  92128
  26.   (619) 451-3752
  27.  
  28.   *********************************************************/
  29.  
  30. #include   <math.h>        /*needed for floating point and math functions*/
  31. #include   <stdio.h>        /* gives fopen, fclose, printf functions */
  32. #include   <strings.h>        /* string manipulation functions like strcpy and strcat */
  33. #include <EventMgr.h>         /* included for Mac event traps only */
  34. #include <MacTypes.h>         /* included for Mac event traps only */
  35.  
  36. #define        NUMROWS         10         /* number of neurodes = NUMROWS*NUMCOLS */
  37. #define       NUMCOLS         10
  38. #define        PATSIZE         2         /* keep input size small for easy plotting */
  39. #define    NUMITERS      2000     /* total iterations to run */
  40. #define        SAVECOUNT      100          /* when to save net (should be about 0.1*NUMITERS) */
  41. #define        ADJNBORS      500          /* how often to lower neighborhood size (about 0.2*NUMITERS) */ 
  42.  
  43. /*-----------------------------------------------------------------------------------------------------------
  44.    Global storage
  45. -------------------------------------------------------------------------------------------------------------*/
  46. float        inpatt[PATSIZE][NUMITERS]                ; /* input pattern data   */
  47. float        weights[NUMROWS][NUMCOLS][PATSIZE]           ; /* weight storage        */
  48. double      gaindelta = 0.010                                     ; /* gain decrease increment (about 0.1*initial gain)*/
  49. int        winrow,wincol                                        ; /* winner neurode        */
  50. int        iteration                                                ; /* iteration count       */
  51. int        neighborsize = 3                             ; /* initial neighbor size */
  52. double      gain = 0.25                            ; /* initial gain value  */
  53.  
  54. FILE    *fopen(),*savefile                            ;/* output file to save network state */
  55. char    *path =   "Kohonen save"                        ;/* file storage for network state */
  56.  
  57. /*-----------------------------------------------------------------------------------
  58.  
  59.    GET_RAND_INPUT -- sets a new random input pattern into inpatt[][]
  60.    This routine should be modified to force different patterns of
  61.    input data (different probability density functions).  One way
  62.    to do this is to force the random values to fit inside various
  63.    envelopes which are more restrictive than a simple +1..-1 range.
  64.    For example, this version forces the data to fit inside a
  65.    circle of radius 0.5, centered about (0.25,0).  It also is 
  66.    set up so that points in the 1st and 2nd quadrants will be
  67.    more likely than points in the 3rd and 4th quadrants.  Thus
  68.    we have a non-uniform probability density function for the
  69.    network to model.
  70.    It should be noted that this routine is strictly arbitrary and can
  71.    be modified or replaced by any other input pattern generator
  72.    you like.
  73.    No parameters, no return
  74.    
  75. -------------------------------------------------------------------------------------*/
  76. get_rand_input()
  77. {
  78.  double  size;
  79.  double  in0, in1;
  80.  double  factor_1_2, factor_3_4;
  81.  double  quad0, quad1;                /* these are multiplicative factors (+1, -1) to determine the quadrant of the
  82.                              input values.  quad0 multiplies the x-coord, quad1 multiplies the y-coord. */ 
  83.  int          i,iter;                /* iteration number */    
  84.      
  85.      factor_1_2 = 0.6;         /* these are relative likelihoods of being in
  86.                                            the 1st/3rd quadrant vs. the 2nd/4th quadrant */
  87.      factor_3_4 = 0.3;
  88.      for (iter=0; iter<NUMITERS; iter++)
  89.      {
  90.         in0 = ((double) (rand())); 
  91.         in1 = ((double) (rand()));
  92.     
  93.         /* force the input pattern to lie in a non-uniform pdf */
  94.         in0 *= factor_1_2; /* first multiply by weighting factor */
  95.         in1 *= factor_3_4;
  96.     
  97.         /* we 'flip a weighted coin' to decide if this input will be
  98.            in 1st/2nd or 3rd/4th quadrants. */
  99.         if (in0 > in1)
  100.         {
  101.             /* 1st/2nd quadrant won, so we only need to decide which
  102.                quadrant (1 or 2) to put it in */
  103.         
  104.             if (2*in1 < in0) 
  105.             {  /* quadrant 1 */
  106.                 quad0 = 1.0;        /* x positive */
  107.                 quad1 = 1.0;        /* y positive */
  108.             }
  109.             else
  110.             {   /* quadrant 2 */
  111.                 quad0 = -1.0;        /* x negative */
  112.                 quad1 =  1.0;        /* y positive */
  113.             }
  114.         }
  115.         else
  116.             {
  117.              /* 3rd/4th quadrant won, so we need to decide which of
  118.                 quadrant 3 or 4 to put it in */
  119.              if (2*in1 < in0)
  120.              {  /* quadrant 3 */
  121.                  quad0 = -1.0;        /* x negative */
  122.                  quad1 = -1.0;        /* y negative */
  123.              }
  124.              else
  125.              {  /* quadrant 4 */
  126.                 quad0 = 1.0;        /* x positive */
  127.                 quad1 = -1.0;        /* y negative */
  128.             }
  129.         } /* end else quadrant choice */
  130.         
  131.         /* we can now get an input pattern.  The rand() returns a
  132.            random integer 0..32767.  Dividing by 32768 forces a 
  133.            range of 0 to 1.  We multiply by the quadrant factor
  134.            to put it in the correct quadrant. */
  135.     
  136.         inpatt[0][iter] = quad0 * (((float) (rand())) / 32768.0);
  137.         inpatt[1][iter] = quad1 * (((float) (rand())) / 32768.0);
  138.     
  139.         /* normalize the input vector */
  140.         size = (double)(inpatt[0][iter] *inpatt[0][iter] + inpatt[1][iter] * inpatt[1][iter]);
  141.         size = sqrt(size);
  142.         inpatt[0][iter] /= (float)size;    
  143.         inpatt[1][iter] /= (float)size;
  144.         printf(" .");
  145.     }  /* end for each pattern */    
  146.     return;
  147. }
  148.  
  149. /*---------------------------------------------------------------------------------------------------
  150.  
  151.     GET_WINNER  -- sets the winner into winrow,wincol
  152.     curr_iter parameter:  current iteration number
  153.     
  154. -----------------------------------------------------------------------------------------------------*/
  155. get_winner(curr_iter)
  156. int    curr_iter;
  157. {
  158.  int        row,col,pat;
  159.  double    mindist,dist;
  160.  double    diff;
  161.  
  162.      mindist=0.0;
  163.      for (row=0; row<NUMROWS; row++)
  164.      {
  165.          for (col=0; col<NUMCOLS; col++)
  166.          {
  167.              dist = 0.0;
  168.              for (pat=0; pat<PATSIZE; pat++)
  169.              {    
  170.                  /* compute distance of this neurode's weight vector to input pattern 
  171.                      we don't care about taking the square root here -- we only want relative
  172.                      distance */
  173.                  diff = (inpatt[pat][curr_iter]) - weights[row][col][pat];
  174.                  dist += diff*diff;
  175.              }
  176.              
  177.              if ( (row==0) && (col==0) )    /* first neurode is always the winner so far! */
  178.              {
  179.                  mindist = dist;
  180.                  winrow  = row;
  181.                  wincol  = col;
  182.              }
  183.              else
  184.              {
  185.                  if (dist < mindist)     /* there's a new winner 
  186.                                       notice that in a tie, the first neurode found 
  187.                                       will be the winner always */
  188.                  {
  189.                      winrow = row;
  190.                      wincol = col;
  191.                      mindist = dist;
  192.                  }  /* end if new winner */
  193.                  
  194.              }  /* end else not first neurode */
  195.          }  /* end for col */
  196.      }  /* end for row */
  197.      return;
  198.  } /* end get_winner */
  199.  
  200. /*----------------------------------------------------------------------------------------------------------
  201.  
  202.     ADJUST_WTS -- adjust the weights of the neurodes in the winning
  203.                       neighborhood
  204.     curr_iter parameter:  current iteration parameter
  205.     
  206. ------------------------------------------------------------------------------------------------------------*/
  207.  
  208. adjust_wts(curr_iter)
  209. int        curr_iter    ;  /* current iteration number */
  210. {
  211.  int    minrow,mincol,maxrow,maxcol;
  212.  int    row,col;
  213.  
  214.  /* neighborhood of winner is determined by 'neighborsize' 
  215.     and is allowed to go to 0 -- i.e., only the winner is 
  216.     adjusted in that case.*/
  217.     
  218.      minrow = winrow - neighborsize;
  219.      maxrow = winrow + neighborsize;
  220.      mincol = wincol - neighborsize;
  221.      maxcol = wincol + neighborsize;
  222.  
  223.      /* make sure that the neighborhood is constrained to be
  224.         inside the network! */
  225.         
  226.      if (minrow < 0)
  227.          minrow = 0;
  228.      if (maxrow >= NUMCOLS)
  229.          maxrow = NUMCOLS - 1;
  230.      if (mincol < 0)
  231.          mincol = 0;
  232.      if (maxcol >= NUMCOLS)
  233.          maxcol = NUMCOLS - 1;
  234.  
  235.     /* we have a good neighborhood, so adjust the weights */
  236.     
  237.      for (row=minrow; row<= maxrow; row++)
  238.      {
  239.          for (col=mincol; col<= maxcol; col++)
  240.          {
  241.              weights[row][col][0] += 
  242.                      gain*(inpatt[0][curr_iter]-weights[row][col][0]);
  243.              weights[row][col][1] += 
  244.                      gain*(inpatt[1][curr_iter]-weights[row][col][1]);
  245.          } /* end for col */
  246.      } /* end for row */
  247.      
  248.      return;
  249.  } /* end adjust_wts */
  250.  
  251.  /*------------------------------------------------------------------------------------------------
  252.  
  253.      DISP_ITERATION -- show iteration number on screen
  254.      Parameter = iteration number
  255.      No return
  256.  --------------------------------------------------------------------------------------------------*/
  257.  disp_iteration (itnum)
  258.  int    itnum;
  259.  {
  260.      printf("\n  Iteration count:  %4d",itnum);
  261.      return;
  262.  } /* end disp_iteration */
  263.  
  264.  
  265.  
  266.  /*-------------------------------------------------------------------------------------------------
  267.  
  268.      INITIALIZE -- initialize weights and random number generator
  269.      No parameters, no return
  270.      
  271.  ---------------------------------------------------------------------------------------------------*/
  272.  initialize()
  273.  {
  274.   unsigned int    seed;
  275.   int            row, col, patt;
  276.   double        number;
  277.   float        sum;
  278.   int            i;
  279.   
  280.   /* initialize random number generator */
  281.   printf ("\n Please enter a random number seed (unsigned integer):");
  282.   scanf ("%d",&seed);
  283.   srand(seed);
  284.   /* initialize weight matrix with random values */
  285.   printf ("\n Initializing weight matrix ");
  286.   for (row=0; row<NUMROWS; row++)
  287.   {
  288.       for (col=0; col<NUMCOLS; col++)
  289.       {    
  290.           sum = 0.0;
  291.           for (patt=0; patt < PATSIZE; patt++)
  292.           {
  293.               number = (double) rand();
  294.               weights[row][col][patt] = (float)number / 32768.0;
  295.             sum += (weights[row][col][patt]*weights[row][col][patt]);
  296.           } /* for each pattern element */
  297.           /* now normalize the weights for this neurode */
  298.           sum = (float)(sqrt((double)sum)); 
  299.           for (patt=0; patt<PATSIZE; patt++)
  300.           {
  301.               weights[row][col][patt] /= (float)sum;
  302.           }
  303.       } /* end for each column */
  304.       printf(" .");                 /* marker on screen to let operator know you're running */
  305.   } /* end for each row */
  306.   return;
  307. } /* end initialize */
  308.  
  309. /*--------------------------------------------------------------------------------------------------
  310.  
  311.     SAVE_DATA(iteration)
  312.     This routine will dump the current state of the weight 
  313.     vectors to a file where they can later be plotted.
  314.     The format of the data file is:
  315.     
  316.     Iteration number
  317.     neurode 0,0 weights 
  318.     neurode 0,1 weights
  319.     neurode 0,2 weights
  320.     ...
  321.     neurode 1,0 weights
  322.     ...
  323.     neurode 9,9 weights
  324.     Iteration number
  325.     ...
  326.     
  327.     The routine or program reading this file must be told the
  328.     number of neurodes in each row and column of the network,
  329.     and the number of elements in the input pattern.
  330.     This routine is called every SAVECOUNT iterations,
  331.     which should be set so there will be about 10 data dumps
  332.     every time this program is run.
  333.     By looking at the data dumps, we will be able to see the
  334.     self organization process.
  335.     
  336.     Parameter iteration is the current iteration count for this
  337.     save.
  338. -----------------------------------------------------------------------------------------------------*/
  339. save_data(count)
  340. int    count;
  341. {
  342.     int        row, col, wt;
  343.     char        *savepath;
  344.     char        *strcpy();
  345.     char        *strcat();
  346.     char        digits;
  347.     char        itoa();
  348.  
  349.     switch (count)
  350.     {
  351.         case 0:  savepath = strcpy(savepath,"Kohonen Iter 0"); break;
  352.         case 100:  savepath = strcpy(savepath,"Kohonen Iter 100"); break;
  353.         case 200:  savepath = strcpy(savepath,"Kohonen Iter 200"); break;
  354.         case 300:  savepath = strcpy(savepath,"Kohonen Iter 300"); break;
  355.         case 400:  savepath = strcpy(savepath,"Kohonen Iter 400"); break;
  356.         case 500:  savepath = strcpy(savepath,"Kohonen Iter 500"); break;
  357.         case 600:  savepath = strcpy(savepath,"Kohonen Iter 600"); break;
  358.         case 700:  savepath = strcpy(savepath,"Kohonen Iter 700"); break;
  359.         case 800:  savepath = strcpy(savepath,"Kohonen Iter 800"); break;
  360.         case 900:  savepath = strcpy(savepath,"Kohonen Iter 900"); break;
  361.         case 1000:  savepath = strcpy(savepath,"Kohonen Iter 1000"); break;
  362.         case 1100:  savepath = strcpy(savepath,"Kohonen Iter 1100"); break;
  363.         case 1200:  savepath = strcpy(savepath,"Kohonen Iter 1200"); break;
  364.         case 1300:  savepath = strcpy(savepath,"Kohonen Iter 1300"); break;
  365.         case 1400:  savepath = strcpy(savepath,"Kohonen Iter 1400"); break;
  366.         case 1500:  savepath = strcpy(savepath,"Kohonen Iter 1500"); break;
  367.         case 1600:  savepath = strcpy(savepath,"Kohonen Iter 1600"); break;
  368.         case 1700:  savepath = strcpy(savepath,"Kohonen Iter 1700"); break;
  369.         case 1800:  savepath = strcpy(savepath,"Kohonen Iter 1800"); break;
  370.         case 1900:  savepath = strcpy(savepath,"Kohonen Iter 1900"); break;
  371.         case 2000:  savepath = strcpy(savepath,"Kohonen Iter 2000"); break;
  372.  
  373.         default:    savepath=strcpy(savepath,"Kohonen Iter");break;
  374.     }
  375.     savefile = fopen(savepath,"w");         /* open an output text file */
  376.     
  377.     fprintf(savefile,"\n Iteration count %d",count);   /* iteration count */
  378.     for (row=0; row<NUMROWS; row++)
  379.     {
  380.         for (col=0; col<NUMCOLS; col++)
  381.         {
  382.             fprintf(savefile,"\n"); /* new line in file */
  383.             for (wt=0; wt<PATSIZE; wt++)
  384.             {
  385.                 fprintf(savefile,"%8.3f",weights[row][col][wt]);
  386.             } /* end for each weight in neurode */
  387.         }  /* end for each neurode in column */
  388.     }  /* end for each neurode in row */
  389.     fclose(savefile);
  390.     return;
  391. }
  392.  
  393.  /*----------------------------------------------------------------------------------------------------
  394.  
  395.      Main program starts here
  396.      
  397.  ------------------------------------------------------------------------------------------------------*/
  398.  
  399.  main()
  400.  {
  401.   int    i;
  402.  
  403.   /* The next two definitions are Mac unique defined in Mac Toolbox 
  404.       You should substitute whatever code is necessary to allow your computer
  405.       to interrupt processing on operator command.  This might be a break or
  406.       other operating system command which you can issue at any time during
  407.       execution of the program.  The idea is to allow you to "bail out" at will */
  408.       
  409.   int            mouse_pressed;         /* used to check for mouse press */
  410.   EventRecord  *event;               /* used to check for mouse press */
  411.     
  412.   /* Set event mask to screen for either mouse pressed or released */
  413.  
  414.   mouse_pressed = BitOr (mouseDown, mouseUp);      /* Macintosh unique code */
  415.  
  416.      initialize();
  417.       get_rand_input();
  418.       save_data(0);                         /* save the initial weights */
  419.       printf("\n  Initial neighborhood size is %d",neighborsize);
  420.       printf("\n  Initial gain is %8.3f",gain);
  421.       for (i=1; i<=NUMITERS; i++)
  422.      {
  423.          iteration = i;
  424.          disp_iteration(iteration);
  425.          get_winner(iteration);
  426.          adjust_wts(iteration);
  427.          if ((iteration % SAVECOUNT) == 0)
  428.          {    
  429.              printf("\n Saving current network state for iteration %d...",iteration);
  430.                  save_data(iteration);            /* dump wts to file for later plotting */
  431.              gain -= gaindelta;               /* lower the gain slightly */
  432.              printf("\n Gain now equals %8.3f",gain);
  433.             /* having the neighborhood size adjustment inside the "data save count" loop
  434.                only works if ADJNBORS is an integral number times SAVECOUNT!  */
  435.  
  436.              if ((iteration % ADJNBORS) == 0)
  437.              {
  438.                  neighborsize -= 1;
  439.                  if (neighborsize < 0) neighborsize = 0;
  440.                  printf("\n Neighborhood size now equals %d",neighborsize);
  441.              } /* end every ADJNBORS iterations */
  442.          } /* end every SAVECOUNT iterations */
  443.  
  444.          /* The following is Macintosh unique code to allow an easy bailout */
  445.         /* Any mouse click will stop the program.  We just force the loop variable
  446.             to be too large to continue the loop. */
  447.         
  448.         if (GetNextEvent(mouse_pressed,event)==TRUE)
  449.             i=NUMITERS+1;
  450.  
  451.      } /* end for NUMITERS iterations */
  452.      printf("\n\n Training is complete.  Saving last state of network now...");
  453.      save_data(iteration); /* dump the last state of the network */
  454.      return;
  455.  } /* end main */